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Abstract 

In  this  report  we  examine  the  Bayesian  method  for  testing  for  compliance 
to  a  given  threshold  studied  by  Nicholson,  Mensing  and  Gray.  It  is  noted  that 
although  this  test  and  accompanying  confidence  intervals  are  valid  for  a  single 
»  event,  it  is  incorrect  to  apply  it  or  the  confidence  intervals  to  repeated  events 

at  the  same  site  unless  the  number  of  calibration  events  is  large.  Since  in  any 
foreseeable  future  the  number  of  calibration  events  is  likely  to  be  small,  this 
report  studies  the  applicability  of  the  Bayesian  test  in  this  case.  The  results 
suggest  that  in  many  instances  the  Bayesian  method  examined  here  should  be 
used  on  repeated  events  with  caution  if  the  number  of  calibration  events  is  less 
than  three. 
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1  Introduction 


Over  the  last  few  years  much  of  the  interest  in  yield  estimation  and  threshold 
test  ban  treaty  monitoring  has  shifted  to  the  problem  of  properly  monitoring 
yields  that  are  somewhat  smaller  than  the  current  test  ban  limit  of  150  Kt. 
As  a  result  of  this  interest  in  smaller  yields  it  has  become  more  important  to 
include  the  effects  of  unknown  slope  (in  the  standard  magnitude/yield  relation) 
on  estimated  yields,  associated  confidence  intervals,  and  related  hypothesis  tests. 

The  most  popular  approach  for  addressing  this  problem  thus  far  has  been 
through  the  Baysian  methodology.  See  W.  L.  Nicholson,  R.  W.  Mensing  and  H. 
L.  Gray,  or  R.  H.  Shumway  and  Z.  A.  Der  for  example.  In  each  of  these  papers 
the  authors  make  use  of  prior  distributions  on  the  parameter  spaces  to  obtain 
estimates  of  yield,  confidence  intervals  for  yield,  threshold  type  test  of  hypotheses, 
and  associated  F-numbers  which  allow  for  errors  in  estimating  geological  bias  and 
slope  as  well  as  several  other  unknown  parameters.  Although  such  results  are 
exactly  what  was  needed  in  one  sense  they  present  a  problem  in  another.  That 
is,  although  the  confidence  intervals  and  hypothesis  tests  are  valid  when  related 
to  a  single  event  from  all  possible  parameter  configurations  they  do  not  represent 
such  intervals  or  hypothesis  tests  when  applied  repeatly  to  a  fixed  test  site  (This 
will  be  explained  in  detail  in  section  4).  This  problem  was  noted  by  Fisk,  Gray, 
McCartor  and  Wilson  (1991)  for  the  case  where  the  slope  is  known. 

In  this  report  we  examine  the  current  Baysian  approach  to  yield  estimation 
from  several  practical  aspects.  That  is  we  consider: 

1.  The  power  of  the  tests  for  several  different  parameter  configurations  and 
yield  training  sets. 

2.  The  maximum  benefit  of  previous  no  yield  data  regarding  its  contribu¬ 
tion  to  increasing  the  power  or  decreasing  the  F-number. 

3.  The  actual  error  rate  or  confidence  interval  (Cl)  that  results  when  these 
Baysian  tests  or  Cl’s  are  applied  to  repeated  tests  at  the  same  site. 

Item  3  is  of  special  interest  if  the  number  of  calibration  events  is  small  and  the 
particular  test  site  is  an  anomaly,  i.e.  a  site  whose  parameters  differ  substantial 
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from  their  corresponding  Bayesian  means.  We  shall  refer  to  our  investigation  of 
item  3  as  a  robustness  study. 

2  Notation  and  Background 

Let  Yj  denote  the  jth  yield  at  a  given  test  site  and  let  my  denote  the  ith 
magnitude  associated  with  the  jth  yield, 

™ij  =  Ai  +  BiW0j  +  ey  ( 1 ) 

*  =  1,2 ,*-*,p  and  j  =  1,2, •••,«,  where  W(y  =  logl^-  -  logFo  =  Wj  - 
Wq ,  with  Wo  given  and  the  ey  represent  random  errors.  Further  let  A  = 
(■^1,* **> Ap)\  B  =  and  ej  =  (ey,e2j, •  •  • , epj)'  where  the  prime 

denotes  transpose,  and  the  ej  are  normal  random  vectors  with  mean  (0, 0,  •  •  • ,  0/ 
and  known  variance  Ee.  We  can  now  write  (1)  in  the  matrix  form 

mj  =  A  +  BWoy  +  ey  (2) 

In  the  model  defined  by  Equation  (1)  A  and  B  are  vectors  of  parameters 
that  depend  on  the  test  site  and  the  particular  magnitude  being  considered.  For 
example  my  may  refer  to  the  jth  mj  value  while  m^j  might  be  the  jfth  m^ 
value.  Ideally  A  and  B  in  (2)  would  be  known.  This  is  in  general  not  the  case. 

However  there  may  be  sufficient  information  regarding  A  and  B  to  restrict 
their  possible  values.  That  is,  it  is  arguable  that  one  can  reasonably  impose 
a  probability  distribution  on  A  and  B  a  priori.  This  is  in  fact  the  reasoning 
that  leads  to  a  Bayesian  approach  to  the  problem.  Specifically  we  suppose  that 
P  =  (A',  B')'  has  a  prior  normal  distribution  with  known  mean  and  covariance 
Up.  In  the  future  we  will  denote  this  by  (3  rsj  Therefore  in  Equation 

(2)  we  no  longer  treat  A  and  B  as  fixed  parameters  but  as  random  variables  or,  if 
you  like,  “parameters”  which  take  on  their  possible  values  with  some  probability. 

Now  suppose  n  calibration  events  are  available,  i.e.  n  events  at  a  given  site 
for  which  the  yields  Wj  are  known  (or  at  least  known  sufficiently  well  that  we 
can  neglect  the  errors  in  the  observed  Wj).  Then  we  can  determine  a  compliance 
test  and  its  associated  F-number  which  properly  integrates  the  information  in  the 
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prior  distribution  with  the  data  from  the  calibration  events.  This  is  the  subject 
of  the  next  section. 

3  A  Bayesian  Test  of  Compliance 


In  order  to  determine  a  test  for  compliance  which  makes  use  of  prior  informa¬ 
tion  regarding  / 3  and  the  calibration  events,  we  need  to  determine  the  probability 
density  function  (pdf)  for  m  =  mn+i  given  mi,  m2,  •  •  • ,  m„.  We  will  denote  this 
pdf  by  /(m|  it5n),  where  itfn  =  (mi  ,  m2,  •  •  • ,  mn)'.  Given  n  events  for  which  the 
yields  are  known  we  wish  to  develop  a  compliance  test  for  an  (n  +  l)st  event  for 
which  the  yield  is  unknown. 

Note  that  the  model  in  (2)  can  be  written  in  the  form 


where 


D  j  = 


nij  = 

4"  j 

=  1,2, •••,«, 

/I 

0  •  •  •  0 

W0j  0 

0  \ 

0 

1  •••  0 

•  .  t  • 

0  w0j 

0 

\0 

0  •••  1 

0  0 

*  •  • 

•  •  •  Wq  j  / 

(3) 


—  (1,  Woj)  ®  I p 


and  <g)  denotes  the  kronecker  product. 

Case  1:  0  known 

The  problem  is  a  simple  one  when  0  is  known  since  in  that  event  m  is 
independent  of  the  previous  mi,m2,  •  •  •  ,mn,  i.e.  Up  =  0.  The  hypothesis  Hq, 
to  be  tested  is 


Ho:  W  <  WT 

against  (4) 

Hi  :  W  >  WT , 

where  W  =  Wn+i.  If  we  shorten  the  notation  T^o.n+l  Wc>  *-e-  take  Wo)n+l  = 
Wc  we  can  write  the  hypothesis  test  in  Equation  (4)  in  the  form 

H0  :  Wc<  WcT 


Hi  :  Wc  >  WcT , 
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where  Wcp  =  Wp  —  Wq.  In  this  case  /(m|mn)  =  /( m).  Now  let 

P 

=  r«m»’  (5) 

*=1 

where  the  rt-  are  known  weights  with  0  <  r,-  <  1  and  £?=1  r*  =  1.  It  is  well  known 
that  if  ej  ~  N( 0, Ee),  then  m  ~  JV(D0,  Ee)  where  D  =  (1,  Wc)®Ip,  and  it  then 
follows  at  once  that  mr  ~  N{v>'D0, r'Eer),  where  r  =  (r1,r2,---, rp)'.  Therefore, 
xmder  Hq,  we  take  Wc  =  Wcp  so  that  a  test  of  the  hypothesis  in  Equation  (4)  at 
the  lOOc*  percent  significance  level  is  given  by  the  following  rule 

Reject  Ho  if  mr  >  Tia,  (6) 

where 

Tia  =  r'DT0  +  zQ  vVEer,  ( 7) 

Dr  =  (1>  WcT)  0  Ip, 

and  za  is  the  100(1  —  a)th  percentile  point  of  N{ 0, 1)  distribution.  We  shall  refer 
to  the  test  defined  by  the  rule  given  by  Equation  (6)  as  Test  1. 

Case  2:  fi  unknown 

Of  course  0  is  not  known  and  therefore  Test  1  cannot  be  used  in  practice. 
It  does  however  furnish  us  a  base  line  for  comparison  purposes.  What  can  be 
reasonably  assumed,  as  we  have  already  mentioned,  is  that  0  where 

Up  and  Up  are  known.  In  this  case  m  and  nin  are  not  independent  and  therefore 
the  problem  is  a  bit  more  difficult.  It  can  however  be  solved  by  making  use  of 
the  following  theorem,  the  proof  of  which  we  include  in  the  Appendix. 

Theorem  1.  Let  m  =  be  a  p-dimensional  magnitude  related  to 

Wb,n+1  =  Wc  by  the  model  of  Equation  (3).  For  k  =  1,2,  •  •  •  ,n  +  1,  let  = 
Hj=l  m/A,  w»^(jfc)  =  E*=l  Wojmj/k,  and  Wok  =  Ej=l  Woj/k.  Suppose  0 
has  the  prior  density  N(pp,  Ep).  Then  the  probability  density  of  m  given  mn, 
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/(m|n£n),  is  N(n,E),  where 

i-l 


E  =  Ee[Ee-H]  Be, 

/i  =  e{(1,  Wc)  ®  +  {e„+i  ®  (Ee/(n  +  1))}] 


-1 


{e„+i  0  (E./(n  +  1))  }  EjV,J  +  n  (l2  ®  Be ')  (  ^  ) 


and 


H  =  {(1,  We)  ®  Ip}*)/?  [S/9  +  {e„+i  ®  (Ee/(n  +  1))  }] 
.{E„+1®(Be/(n  +  l))}{(^) 

E  -(  1  W°'"+1  V1 

En+1 "  Uo,n+i  E”i/ <•/(«  +  !)  i 


-1 


(8) 


(9) 


(10) 


Given  Theorem  1,  the  problem  is  once  again  trivial  and  we  can  again  write 
down  a  100o%  significance  level  test.  If  mr  is  defined  by  Equation  (5),  then  it 
follows  from  Theorem  1  that  the  pdf  of  mr  given  mn  is  Ar(r'/i,r'Er),  where  (i 
and  E  are  given  by  (8)  and  (9)  respectively.  The  desired  test  is  then 

Reject  Hq  if  mr  >  l2a, 


where 

?2a  =  *'/*  +  *aV r'Er,  (11) 

H  and  E  axe  defined  in  Theorem  1  with  Wc  =  Wct,  and  zQ  is  the  100(1  —  a) 
percentile  point  of  iV(0, 1).  We  shall  refer  to  the  test  of  Equation  (11)  as  Test  2. 

From  Theorem  1  one  can  also  obtain  a  confidence  interval  for  W .  For  a 
general  treatment  of  the  confidence  interval  problem  with  a  Bayesian  prior  see 
R.  H.  Shumway  and  Z.  A.  Der. 


6 


4  A  Constrained  Bayesian  Test 


In  a  recent  report  Nicholson,  Mensing  and  Gray  (1991)  show  how  previous 
magnitude  data  can  be  used  to  define  a  Bayesian  prior  for  even  though  the 
associated  yields  are  not  available.  We  shall  refer  to  such  data  as  “no  yield” 
data  as  before,  and  we  also  assume  that  n  calibration  events  are  available.  In 
this  section  we  consider  the  question,  “What  is  the  maximum  information  that 
can  be  gained  by  this  approach  ?”  In  order  to  accomplish  this  we  will  consider 
the  problem  of  the  previous  section  but  we  let  the  number  of  no-yield  events 
go  to  infinity.  That  is,  we  consider  the  case  where  the  “no  yield”  data  set  is 
sufficiently  large  that  the  parameters  that  are  estimable  from  that  data  can  be 
estimated  without  error,  *.e.  they  are  known.  By  developing  a  test  for  this 
case  and  comparing  its  power  to  Test  2  we  are  able  to  determine  the  maximum 
improvement  in  power  (or  reduction  in  F-number)  obtainable  in  the  approach  of 
Case  2. 

Specifically  we  note  that  the  parameters 

Cj_ i  =  i  =  2,3, •••  ,p 

and  (12) 

hi— 1  = 

do  not  depend  on  yield  and  hence  consistent  estimates  for  c,-_i  and  //;_i  can  be 
obtained  from  the  “no  yield”  data.  Thus  in  this  case  we  take  ct-_i  and  /z,-_ \  as 
known,  t  =  2,  •  •  •  ,p.  Moreover  under  these  constraints  the  model  of  Equation  (3) 
becomes 


mj  =  Pi  +  &LjPi  +  ej,  (13) 

where  fa  =  (Ai, #i)\  hL  =  (0, hU l*2> •**, hp-l)'  and 

D Lj=(  1  C1  C2  “*  Cp_1  \  . 

\W0j  ciW0j  C2W0j  •••  Cp^iWojJ 

Thus  the  original  2 p  dimensional  parameter  space  for  ft  is  reduced  to  the  2 

dimensional  one  due  to  the  constraints  in  Equation  (12). 


7 


To  determine  a  test  for  the  model  of  Equation  (13)  we  need  the  following 
theorem,  the  proof  of  which  is  included  in  the  Appendix. 

Theorem  2.  Let  m  =  mn+i  be  a  p-dimensional  magnitude  related  to 
Wjj  n+1  =  Wc  by  the  model  of  Equation  (13).  Suppose  has  the 

prior  pdf  JV(jii,Ei).  Then  the  pdf  of  m  given  mn  is  7V(/ic,Ec),  where 


Ec  = 


Ip  ~  ^I,n+lQn+l^i,n+l^e  Ee, 


-1 


flc  =  ML  +  £cEe  ^L,n+ lQ„+lR"’ 


-1 


and 


n+1 

Qn+1  =  +  y>^e  ^ij), 

3= 1 

R„  =  E^Vl  +  DIiEe  !(mi  -  VL), 

3= 1 


D 


ci  c2 


cp— 1 


L,n+1 


Wc  ciWc  C2WC  •••  Cp-xWc. 


From  Theorem  2  it  follows  that  mr  ~  N(r'pc,T'%cr)  and  hence  a  100(o)% 
significance  test  of  the  hypothesis  in  Equation  (4)  is  given  by  the  following  rule: 

Reject  Hq  if  mr  >  T^Ql 

where 


T3a  =  r'nc  +  za  y/r%*,  (14) 

and  za  is  the  100(1  —  a)  percentile  point  of  a  N( 0, 1).  We  shall  refer  to  the  test 
of  Equation  (14)  as  Test  3. 

5  Power  Curve  Comparisons 

In  order  to  assess  the  impact  of  imposing  the  prior  information,  we  compare 
the  power  of  the  following  tests: 
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Test  1  :  a  test  of  hypothesis  based  on  the  assumption  that  the  population 

parameters  are  known. 

Test  2  :  a  test  of  hypothesis  based  on  the  unconstrained  Bayesian  ap¬ 

proach  and  the  assumption  that  the  parameters  are  unknown. 

Test  3  ’■  a  test  of  hypothesis  based  on  the  constrained  Bayesian  approach 
and  the  assumption  that  the  population  parameters  are  un¬ 
known. 

The  power  at  W  is  given  by 


Power  (W)  =  P(mr  >  Ta|mn,  W). 

Also  the  F-number  of  the  test  is  given  by 

p  _  10WF-WT , 

where  WF  is  the  value  of  the  log  yield  at  which  the  power  is  0.5. 

Since  we  specified  the  critical  values  of  Test  1,  Test  2,  and  Test  3  in  (7),  (11), 
and  (14),  respectively,  it  is  easy  to  show  that  the  power  of  Tests  1,  2,  and  3  are 

Power(W)i  =  1  -  $  (ztt  +  ,  (15) 

Power(W)2  =  1  -  $  ( —  ~  ,  (16) 

\  VrSWr  / 

Power(W)3  =  1  -  $  ( ~  Pew)  +  ^avVSc r 
\  Vr'^cWr 

where  D  =  (1,  PF)  ®  Ip,  $  is  the  cumulative  distribution  function  of  7V(0, 1),  and 

{PW,  2^},  {Hew,  2cW/}  are  defined  as  in  Theorem  1,  and  Theorem  2,  respec¬ 
tively,  with  W  given. 

From  (16)  and  (17)  it  is  clear  that  the  power  of  Test  2  and  Test  3  depends  on 
the  value  of  mn.  Therefore  in  order  to  compare  with  Test  1  we  generate  two  equiv¬ 
alent  data  sets  for  Test  2  and  Test  3  with  fixed  values  of  {'2e,fi0,'£0,ci,ni,n} 
when  p  =  2.  With  the  known  parameter  and  the  generated  data  sets,  we  com¬ 
puted  the  power  of  Test  1,  Test  2,  and  Test  3  on  the  100  equally  spaced  grid 
values  between  log  150  and  log  300  for  W  from  (15),  (16),  and  (17),  respectively 
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and  W0  =  log  125.  We  ran  this  simulation  20  times  to  get  the  mean  of  the  powers 
for  Test  2  and  Test  3.  Power(W)i,  mean  Power  (W)2,  and  mean  Power(W)3  are 
plotted  on  Figure  1  through  Figure  8  for  various  values  of  {Ee,  Up, %0, c\,  ”}• 

Now  we  summarize  some  findings  from  the  simulation.  As  we  can  see  in  Fig¬ 
ure  1  through  Figure  3,  mean  Power  (W)2  and  mean  Power(W)3  rapidly  converge 
to  Power(W)i  as  n  gets  large.  Similarly  average  F2  and  average  F3  converge  to 
Fi  as  n  grows,  where  Fi,  F2  and  F3  are  the  F-numbers  of  Test  1,  Test  2,  and 
Test  3,  respectively. 

The  relatively  better  performance  of  Test  3  over  Test  2  is  observed  regardless 
of  the  values  of  c\  in  Figure  1  and  Figure  4.  However,  Figure  5  and  Figure  6 
show  that  the  overperformance  of  Test  3  against  Test  2  diminishes  as  the  standard 
deviations  of  A2  and  B2  (°rA2, &b2)  decrease  to  those  of  A\  and  B\ ,  respectively. 
Figure  7  and  Figure  8  show  the  same  phenomenon  as  ae2  becomes  small  enough 
to  be  similar  to  atl.  Thus  it  would  appear  that  if  the  values  used  here  for  up,  Up 
and  Ee  are  representative,  additional  no  yield  data  would  be  of  little  value. 

6  Robustness 

In  the  previous  sections  we  have  developed  a  test  of  the  hypothesis  of  com¬ 
pliance  of  the  (n  +  l)st  event  given  n  calibration  events  when  fi  is  unknown.  We 
referred  to  this  as  Test  2.  In  making  use  of  this  test  it  is  important  to  under¬ 
stand  the  nature  of  the  false  alarm  rate  or  significance  level  a.  Possibly  the  best 
way  to  interpret  a  is  to  think  through  a  simulation  for  estimating  a.  In  order 
to  simulate  the  process  one  would  first  generate  fi  from  N(fip,  E^)  and  then, 
given  fi  and  W,-,  t  =  1, ••*,«,  generate  6i,62,***,6n  to  obtain  mn.  Now  let¬ 
ting  Wn+i  =  log  150  and  generating  en+i  to  obtain  mn+i,  one  would  apply  the 
test  and  note  the  decision.  This  simulates  the  senario  of  obtaining  n  calibration 
events  and  one  additional  event  of  unknown  yield.  This  entire  process  would  be 
repeated  a  large  number  of  times  and  the  proportion  of  incorrect  decisions  would 
approach  a.  Table  1  below  describes  the  method.  The  fii  denote  the  values  of 
fi  generated  on  simulation  #  i.  Let  mrtn+\(i)  =  where  m^n+i  is  the 
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(n  +  l)st  magnitude  vector  generated  in  the  ith  simulation,  i  =  1,2,***,/. 


Table  1.  Simulation  procedure  for  estimating  the  false  alarm  rate 


Simulation  #  1 

Simulation  #  2 

Simulation  #  l 

Given 

02,Wl 

-  ,Wn 

Generate 

n»ll,  **•»  mln 

“21,  •••,*“ 2„ 

m/i,  *»ra/n 

Generate 

rai(B+i,  W  =  log  150  m2)„+i,  W  =  log  150 

•  •  •  m/n+i ,  W  =  log  150 

Decision 

Reject  if 

Reject  if 

Reject  if 

mr,n+l(l)  >  ^2a(l) 

™r,n+l(2)  >  T2a(2) 

mr,n+l(0  >  T2a(l) 

Now,  if  we  define  a  random  variable  X  such  that  X  —  1  if  mr  n+i(j)  >  T2a(j), 
and  otherwise  X  =  0,  it  follows  that 

I 

a  =  lim  V''  Xj/l  =  lim  Xi  a.s. 

/-oof-'  3  l-+oo  1 
J=1 

We  should  note,  however,  that  in  practice  the  application  of  these  tests  will 
be  to  events  *  *  *  ,mn+s  at  the  same  site.  That  is,  what  is  needed 

is  essentially  a  test  so  that  the  empirical  false  alarm  rate  or  significance  level 
approaches  a  as  s  gets  large  rather  than  as  l  gets  large.  We  shall  refer  to  the 
sample  false  alarm  rate  as  s  — ►  oo  as  the  “actual  significance  level”  or  “actual 
false  alarm  rate”  and  denote  it  by  a(mn|n,/?).  Thus 

^(mnln,/?)  —  Pip^n+s  >  ^2al Wn+i  =  log  150,^)  (IS) 

=  Km  (#  of  mn+k  >  T2a)/k. 

k— KX> 

It  can  be  shown  that 

lim  a(m„| n,B)  =  a. 

Thus  when  n  is  large,  a(mn|n, /))  «  a  regardless  of  the  observed  value  of  fi. 
However  in  most  instances  n  will  be  small  and  therefore  the  question  which 
arrises  is,  “How  robust  is  a  to  small  values  of  n  and  unusual  values  of  ft,  i.e. 
values  of  far  removed  from  /i^j  ?”  That  is,  “How  close  is  a  to  a(mn|n, fi), 
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the  actual  false  alarm  rate,  when  n  is  small  and  is  substantially  different  from 
P/J?” 

In  addition  to  the  “actual  false  alarm  rate”  we  need  to  obtain  the  probability 
of  rejecting  Hq  as  s  — i ►  oo.  We  shall  refer  to  this  as  the  “actual  power”  or  the 
“actual  probability  of  detection”,  and  denote  it  by  P(W\n,f3).  Thus 

P(W\nJ)  =  P(mn+S  >  T2a\mnJ,Wn+8  =  W)  (19) 

=  Hm  (#  of  mn+k  >  T2a)/k. 

k-*oo 

Then  it  also  can  be  shown  that 


lim^PiWlnJ)  =  Power(W). 

From  (18)  and  (19)  it  is  clear  that  a(mn|n,/?)  and  P(W\n,fi)  depend  on  inR. 
Thus  for  every  sample  of  rnn  these  quantities  will  be  different.  We  can  however 
estimate  F?[a(mn|n,)9)]  and  E[P{W\n,fi)]  for  various  values  of  fl  and  n.  This  is 
the  topic  of  the  remaining  portion  of  this  section. 

In  order  to  investigate  the  robustness  of  the  actual  false  alarm  rate, 
a(n?n|n,  /?),  a  small  simulation  was  performed  for  a  variety  of  values  of  n  and  fi. 
Specifically,  taking  p  =  2,  /i^  =  =  (4, 4, 1,1)',  aAl  =  aA^  = 

=  °B2  =  0.05,  PA  =  PB  =  0.5,  PAB  =  o,  aei  =  at2  -  0.05,  pe  =  0.5,  W  = 
log  150  and  Wq  =  log  125,  we  considered  the  cases 


P  =  P/}  +  C-(aAl,aA2,0,0)', 


where  (7  =  0,  ±1,  ±2,  for  n  =  1, 2, 3, 5, 10  and  100. 

For  each  case  a  value  of  mn+i  was  obtained  10,000  times  (or  equivalently 

mn-fi,  *  =  1*  •  •  ’  >  10, 000  was  obtained)  and  a(ntn|«,^)  was  estimated  by 

#  of  rejections 


d(Hln|n,£)  = 


10,000 


(20) 


As  already  noted  o(mn|n,  fi)  depends  on  mn  and  clearly  the  same  is  true  about 
&(itiln|rc,  fi).  Therefore  a  reasonable  measure  of  the  robustness  of  Test  2  when  n 
is  small  is  the  I7[a(mn|n,)9)]  =  fiQ.  To  obtain  an  estimate  of  /za,  for  each  case 

we  generated  20  repetitions  of  d(rrln|n,/f),  t.e. 

,  20 


The  results  of  these  simulations  are  given  in  Table  2  for  a  =  0.025.  It  is  worth 
noting  the  relatively  large  standard  deviation  of  a(m„|n,j5).  In  view  of  the  val¬ 
ues  of  fiQ  one  can  conclude  that  the  distribution  of  a(mn|n,/9)  is  quite  skewed 
to  the  right  or  at  least  contains  some  extreme  values  on  the  right  side.  That 
is,  values  of  d(mn|n,  /?)  much  larger  than  jxa  are  more  frequent  than  values  of 
d(m„|n,/?)  less  than  p.Q,  or  substantially  larger  values  of  d(m„|n,/?)  than  p,Q 
may  not  be  unusual.  Since  d(nin|n,  fi)  is  obtained  from  10,000  repetitions,  it  fol¬ 
lows  that  d(mn|n,/?)  «  o(mn|n, 0).  So  similar  remarks  can  be  made  regarding 
»(mn|n,/9).  The  result  of  this  is  that  Table  2  presents  these  results  in  a  conser¬ 
vative  way  since  most  people  would  interpret  the  mean  as  a  typical  value  of  the 
false  alarm  rate.  What  we  are  cautioning  here  is  that,  in  fact  false  alarm  rates 
substantially  larger  than  the  mean  values  shown  in  Table  2  will  be  much  more 
common  than  in  a  symmetric  distribution.  We  probably  should  have  included 
the  median  in  Table  2,  but  that  was  not  calculated. 

It  should  be  noted  that  if  C  <  0,  the  Bayesian  estimator  of  yield  will  un¬ 
derestimate  yield  and  hence  the  true  false  alarm  rate  will  be  too  small  while  if 
C  >  0  the  estimator  will  overestimate  yield  and  hence  the  false  alarm  rate  will 
be  too  large.  From  inspection  of  Table  2,  it  appears  that  if  we  have  only  1  or  2 
calibration  events,  this  effect  can  be  large,  and  hence  in  this  case  the  Bayesian 
significance  level  or  Cl  may  be  seriously  in  error.  On  the  other  hand  if  n  >  5  the 
method  might  be  considered  adequate,  even  though  for  C  <  0  the  false  alarm 
rate  may  still  be  sufficiently  too  small  that  it  could  very  adversely  effect  the 
power,  i.e.  the  chances  of  detecting  a  violation. 

Power  Considerations 

Figures  1-8  compare  the  power  of  Test  1,  Test  2,  and  Test  3  for  various 
parameter  configurations.  As  in  the  case  of  the  false  alarm  rate,  if  these  parame¬ 
ter  values  are  representative,  little  is  to  be  gained  from  additional  no  yield  data. 
Also,  from  the  comparison  of  the  F-numbers  it  does  not  appear  that  a  great  deal 
is  to  be  gained  by  taking  n  >  2.  Unfortunately  these  rather  pleasant  results  do 
not  uniformly  extend  to  the  actual  power. 
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Figure  9  through  Figure  36  compare  the  “actual”  power  of  Test  2  to  the 
power  of  Test  2,  i.e.  they  compare  P(W|n,  fi)  to  P(W).  The  figures  also  compare 
the  F-number  for  Test  2  to  the  “actual”  F-number.  For  n  <  2  it  is  clear  that 
both  the  power  and  the  F-number  are  seriously  effected  if  C  =  ±2  and  the  same 
is  true  for  C  =  ±1  if  n  =  1.  It  should  be  pointed  out  that  the  small  F-numbers 
associated  with  C  <  0  are  a  result  of  very  large  false  alarm  rates  and  should  not 
be  viewed  as  improved  tests. 

Concluding  Remarks 

In  this  report  we  have  investigated  the  robustness  of  the  Bayesian  method 
(referred  to  as  Test  2)  for  testing  compliance  of  an  observed  yield  to  a  threshold. 
Although  the  simulations  reported  here  were  not  exhaustive,  they  were  adequate 
to  demonstrate  that  the  Bayesian  method  for  testing  compliance  is  probably  not 
satisfactory  if  there  are  only  one  or  two  calibration  events.  Moreover  it  is  highly 
desirable  to  have  five  or  more  calibration  events  to  guarantee  good  agreenent 
with  the  stated  significance  level.  Similar  remarks  could  be  made  regarding  the 
corresponding  confidence  intervals. 

The  consequence  of  these  findings  is  that  if  it  is  unlikely  that  several  cali¬ 
bration  events  will  be  available,  Test  2  and  confidence  intervals  associated  with 
Test  2,  the  Bayesian  tests  and  Cl  discussed  by  Nicholson,  Mensing  and  Gray,  and 
those  introduced  by  Shumway  and  Der  should  be  used  with  care.  In  fact  if  the 
number  of  calibration  events  is  less  than  3  it  would  probably  be  wise  to  consider 
a  constrained  likelihood  method  as  an  alternative  to  the  Bayesian  method,  or,  if 
possible,  the  Bayesian  method  should  be  extended  to  include  the  case  of  several 
events  following  the  calibration  events. 
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Table  2.  Estimate  of  Actual  False  Alarm  Rate  E[a(n,f3)],  a  =  0.025 


C  =  —2 

C=  -1 

<7  =  0 

<7  =  1 

(7  =  2 

n  =  0 

A  a 

0.0000 

0.0000 

0.0023 

0.0488 

0.3165 

st.  dev.  Aa 

0.0000 

0.0000 

0.0001 

0.0006 

0.0011 

st.  dev.  d(mn|n,A) 

0.0000 

0.0000 

0.0004 

0.0025 

0.0050 

n  =  1 

Aa 

0.0006 

0.0034 

0.0149 

0.0532 

0.1418* 

st.  dev.  Aa 

0.0003 

0.0013 

0.0044 

0.0112 

0.0205 

st.  dev.  d(mn|n,/?) 

0.0012 

0.0058 

0.0197 

0.0499 

0.0919 

n  =  2 

Aa 

0.0034 

0.0089 

0.0225 

0.0498 

0.0985 

st.  dev.  A  a 

0.0014 

0.0030 

0.0061 

0.0107 

0.0173 

st.  dev.  d(nln|n,A) 

0.0062 

0.0132 

0.0273 

0.0479 

0.0774 

n  =  3 

Aa 

0.0046 

0.0099 

0.0193 

0.0364 

0.0643 

st.  dev.  /ia 

0.0018 

0.0033 

0.0055 

0.0087 

0.0132 

st.  dev.  d(mn|rc,A) 

0.0079 

0.0146 

0.0246 

0.0389 

0.0592 

n  =  5 

Aa 

0.0076 

0.0121 

0.0192 

0.0283 

0.0425 

st.  dev.  fi a 

0.0019 

0.0029 

0.0042 

0.0058 

0.0081 

st.  dev.  d(m„|n,A) 

0.0084 

0.0128 

0.0190 

0.0259 

0.0363 

n  =  10 

Aa 

0.0119 

0.0155 

0.0201 

0.0255 

0.0316 

st.  dev.  Aa 

0.0017 

0.0025 

0.0030 

0.0034 

0.0040 

st.  dev.  d(mn|n,A) 

0.0075 

0.0110 

0.0132 

0.0152 

0.0179 

n  =  100 

Aa 

0.0229 

0.0230 

0.0234 

0.0238 

0.0240 

st.  dev.  A  a 

0.0015 

0.0013 

0.0016 

0.0014 

0.0015 

st.  dev.  d(rnn|n,/9) 

0.0065 

0.0058 

0.0072 

0.0063 

0.0065 
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*  note:  For  symetric  confidence  intervals  a  100(1  —  2 a)%  two  sided  confidence 
interval  corresponds  to  a  one  sided  a— level  significance  test.  For  example,  for 
Test  1  of  size  0.025,  the  corresponding  two  sided  confidence  interval  is  a  95% 
C.I.  This  suggests  that  if  the  “actual”  significance  level  is  0.14,  the  actual  C.I. 
could  be  a  72%  C.I.  That  is,  if  the  site  geological  bias  is  2<r  greater  than  the 
expected  bias,  ha,  then  even  though  the  Bayesian  significance  level  is  0.025  and 
the  Bayesian  C.I.  is  0.95,  the  actual  significance  level  is  estimated  here  as  0.14 
and  one  would  assume  that  the  actual  two  sided  C.I.  is  around  72%. 
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APPENDIX:  PROOFS 


Proof  of  Theorem  1 

For  the  new  observation  m  related  to  Wc,  the  conditional  pdf  of  m  given 
m„,  /(m|ntn)  is  as  follows: 

/(m|mn)  =  J  fi(m,0\mn)dp 

=  J  f2(^\fiMn)h(P\mn)d0 

=  J  f2(m\P)h(ft\r&n)dfi,  (Al) 

where  /i,/2,  and  fy  axe  the  probability  densities.  The  last  equation  is  obtained 
due  to  the  independence  between  m  and  m„  when  0  is  given.  The  conditional 
distribution  of  0  given  mn  may  be  computed  using  Bayes’  law  as  follows: 


MP\*nn) 


hjm^n  \0) 

f  h{0)L{T&n\0)dp' 


(A2) 


where  h  is  the  prior  density  of  the  parameter  vector  0  and  L  is  the  likelihood 
function  for  the  data  mn,  given  values  of  0.  If  we  assume  ej  are  independent 
multivariate  normal,  then 


law = n  * m#)’ 

3= 1 


where 


Hmj\P)  =  (27r)  p|2e|  1/2exp{  -  -  Dj0)' Ee  -  D;^)}. 

Note 

f2(m\P)L(TZn\P)  =  L(tin+l\0,mn+1  =  m,Wn+1  =  Wc) 
since  ej  axe  independent.  Thus  referring  to  (Al)  and  (A2)  leads  to 


/(m|mn)  =  m>  ^n+1  ~  Wc)d0 

s  KP)ui&„\$)ie 


Thus  if  h(0)  is  available,  /  is  completely  determined. 
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Note 


h(0)L( mn\0)  « 

+  -  Djft'Z-Hmj  -  Djp)}] . 

i= 1 

The  exponential  of  (A3)  is  —1/2  times 

/>  V  +  E  D^D#  -  2^>  +  E 


(A3) 


i=i 


;=i 


+  Y  m'2e  Vj. 
i=l 

Let  Won  =  £j=l  W();7n-  Since  D;  =  (*>  Wbj)  ®  !p> 


(A4) 


E  DjSe'Dj  =  E  ((1.  W'0;)'  ®  Ip)Se  1  ((1,  Wo,)  ®  Ip) 

i=i  j= 1 


-tlf1  ^W) 

,tllW  wl)  i 


-1 


where 


fr[\woj  Wqj 
(won  E"=lW§/nj®(  e 

-If1  )"«m) 

\Vw0„  Ei=l  wljln)  V  'J 

—  |En  ®  ^Ee/n^  |  > 

/I  W0n  ’ 

E  =  I 

”~Won  E"=lW02,/n 


-1 


Let  ih(n)  =  YJj= l  mj/n  and  A^(n)  =  £"=i  Wojmj/n.  Then  it  is  easy  to 
verify  that 

Y  mjfSe  1  Dj  =  (™(n) ’  *!?(»))  (l2  ®  (Ee/n)  )  * 

i=i 
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We  can  rewrite  (A4)  as  follows: 

&  [V  +  {E«®  (Ee/n)  }  ]  P  ~  2  1  +  (*”(n)  ’  mW(n)) 

-1  ” 

•  (l2  ®  (Se/n)  \\fi  +  fi'pZp  1H0  +  Y  m'  Se 1 

j= 1 

=  (fi- Z„)'  [e^1  +  {e„  ®  (Ee/n)  }_1]  (fi  -  Z„) 

“  Zn  [E^  +  {En®  (Se/n)  }  ] Z„  +  V/J  +  Y mjEelmj, 

i=i 

where 


Zn  —  +  ^En  ®  ^Ee/n^  |  j  +  |l2  ®  ^Ee/n^ 


-U  /  m(n) 

.  ^^(n)  /  J 


Since  -(1/2)09  -  Zn)'[E^  +  {E„  ®  (Ee/n)}  *](/?  —  Zn)  is  the  exponential 
of  the  multivariate  normal  density  with  mean  Zn  and  variance  [E^1  +  {En  ® 
(Ee/n)}-1]-1,  and  f  exp[-(l/2)(0  -  Zn /[E^1  +  {En  ®  (Ee/n)}-1](/9  -  Zn)]dp 
is  a  constant,  it  can  be  shown  that 

J  h{p)L(r&n \fi)dfi  oc  exp[-  (l/2){  -  Z'n  [E^1  +  {e„  ®  (se/n)  }_1]z„ 

n 

+  +  Y  mjfSe  lmi}]  • 

i=i 


For  the  new  observation  m,  let  /(m|mn)  be  the  conditional  density  of  m 
given  inn  in  the  unconstrained  case.  Then 


\ _ JK$)L(  mn+i|0,m„+i  =  m,  W0)„+i  =  Wc)dfi 

t ( m  mn )  — - ■ - 

JWWAnWdfi 

oc  exp[  -  (1/2)  {  -  Z'n+1  [E^1  +  {E„+i  ®  (Ee/(n  +  1))  }  *]  Zn+1 
+  Z„  Je^1  +  |e0  0  (lle/nj  |  Jz„  +  m,E“1mJ],  (A5) 
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where 


Zn+l  =  Rn|i  Ep1H0+  {l2®  (Ee/(n  +  l)) 

R-n+1  =  +  {^n+1  ®  (Se/(n+ 1))}  > 

=  (  i  «v+i  r1 

"+1  Wo,»+i  E^;n/("+i)/ 


-1-.  /  ni 


E°+1  Wo,n+l  E"il^ 

n+1 

“(n+l)*  53mi/(n  +  1)’ 

i=i 

n+1 

™W(n+ 1)  =  S  W0m/(n  +  !)- 

i=i 

n+l 

Wo, n+1  =  E  *%/("  +  1), 

J=1 

with  mn+i  =  m  and  Wo, n+1  =  Wc. 
Note 


(n+1) 


mW^(n+l) 


*(n+i)  =  (n/(n  +  l))m(n)  +  (l/(n  +  l))m 

%V(n+l)  =  (n/(n  +  l))miy(n)  +  (wc/{n  +-  l))m. 

Then 

f  m(”+1)  )  =  ("/("  +  !))  (  m<n)  )  +  (l/(»  +  l))  (Mem. 

\mW(n+l)/  \mW(n)/  \WC/ 

Therefore  the  exponential  of  (A5)  is  —1/2  times 

~  M'n+1Rn+iMn+i  +  Z'n  [e^1  +  |En  ®  ^Ee/n^|  jzn  -  /i'E-1/! 


®  m. 
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where 


Mn+i  =  Rn h  +  (n/(n  +  1))  {l2  ®  (Ee/(n  +  1))  } 

E  =  Ee[l!e-H]"1Ee, 

fi  =  e{(1,  Wc)  ®  E^E,  [E^  +  {En+1  ®  (E e/(n  +  1))  }] 

*  {En+1  ®  (Ee/(n  +  1))  }  +  n(l2  ®  EeX)  f 


-I'i  /  m/ 


H  =  {(1,  wy  ®  i,}^  [s^  +  {e„+i  ®  (Ee/(n  +  1))  }]  1 

*  (  En  ■  1  ®  (Ee/(rc  +  1 ) J  |  ^  ®  • 

Since  the  first  three  terms  in  (A6)  are  not  function  of  m,  which  are  constants, 
the  theorem  holds. 

Proof  of  Theorem  2 

Note  the  distribution  of  my  given  fi\  is  the  multivariate  normal  with  mean 
PL  +  &LjPl  and  variance  Ee. 

h(Pi)L(rnn\Pi)  «  exp[  -  (1/2)  j(0i  -  /ii/Ef  1(/9i  -  m) 

n 

+  53(mj  -  VL  -  DLjPi)'2el(™ j  -HL-  Diy/?i)}]  .(A7) 
i=  1 

The  exponential  of  (A7)  is  —1/2  times 

(ft  -  Z„)'  [Er1  +  £  -  Z„)  -  Z’„  [Er1  +  £  (D^E-'Dy)] 

i=  i  j= l 

n 

•  Zn  +  ^lEj-Vi  +  5^(my  -  hl)^- l(mj  -  hl), 
j=  l 


where 


Z„  “  E1 1  +  X)  (Di;Ee  ^-Lj)]  [E1  Vl  +  53  Dl;Ee  !(mj  “  A*Z-)]  • 
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Hence 

Jh(0  l)L(t&n\filWl 

oce*p[-(l/2){-  Z'„  [e^1  +  £  (D^-E-^i,)] Z„ 

3= 1 
n 

+  Mi^Vl  +  -  /ix)}l . 

i= i 

For  the  new  observation  m,  let  /c(m|rrln)  be  the  conditional  density  of  m 
given  m„  in  the  constrained  case.  Then 

,  ,__l=*  x  A(0l)£(mn+i|0i,mn+i  =  m,W0,„+i  =  Wc)dp i 
/c(m  mn)  — - t  ”7 - 

/M/WHUAWi 

n+1 

«  exp[  -  (l/2){  -  Z'n+1  [EJ-1  +  £  (D'^Ee’Di^Jz^, 

3= 1 

+  Z'»  [Er'  +  E  (DiiEe 'Dii)]  Z»  +  (">  -  W)'Ee  *(m  -  W)}]  . 

3= 1 


where  Zn_|.i  is  defined  as  Zn  with  mn+i  =  m.  and 

/I  ci  C2 


DI,n+l  = 


We  ciWc  c2Wc 


cp-lWc 


For  k  =  1,2,  •  •  •  ,n  +  1,  let  R*  =  Ex  Vl  +  ]C}=1  1(mi  ”  Ml))  anc*  ^et 

Qt  •).  Then  with  Z„+1  =R„+Di  n+1E-I(m-W), 

it  is  easy  to  show  that  (A8)  is 

exp[  —  (1/2)  j  —  R'nQjjRn  +  Z„QnZn  +  (m  —  Ml/  [^e  1  —  1^I,n+lQn+l 
•  Di,„+12e-‘] (m  -  pL)  ~  2<QJ*  1Di,„+1E-1(m  -  W)}] 
oc  ezp[-(l/2){(m-/lc)'E;r1(m-/tc)}]> 


where 


Sc  =  Ip  —  Dj^n+lQn+i^L.n+l^e  *J  Ee, 

Me  =  Ml  +  £CE“^ 1Dx)n+iQ^|iR-n- 

The  last  result  is  obtained  because  Qn>  Qn+lj^ri)  and  Zn  are  not  function  of  m. 
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